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MULTIGRID  APPROACHES  TO  NON-LINEAR  DIFFUSION  PROBLEMS  ON 

UNSTRUCTURED  MESHES 


DIMITRI  J.  MAVRIPLIS* 

Abstract.  The  efficiency  of  three  multigrid  methods  for  solving  highly  non-linear  diffusion  problems 
on  two-dimensional  unstructured  meshes  is  examined.  The  three  multigrid  methods  differ  mainly  in  the 
manner  in  which  the  non-linearities  of  the  governing  equations  are  handled.  These  comprise  a  non-linear 
full  approximation  storage  (FAS)  multigrid  method  which  is  used  to  solve  the  non-linear  equations  directly, 
a  linear  multigrid  method  which  is  used  to  solve  the  linear  system  arising  from  a  Newton  linearization  of  the 
non-linear  system,  and  a  hybrid  scheme  which  is  based  on  a  non-linear  FAS  multigrid  scheme,  but  employs  a 
linear  solver  on  each  level  as  a  smoother.  Results  indicate  that  all  methods  are  equally  effective  at  converging 
the  non-linear  residual  in  a  given  number  of  grid  sweeps,  but  that  the  linear  solver  is  more  efficient  in  cpu 
time  due  to  the  lower  cost  of  linear  versus  non-linear  grid  sweeps. 

Key  words,  non-linear,  unstructured,  multigrid 

Subject  classification.  Applied  and  Numerical  Mathematics 

1.  Introduction.  The  iterative  solution  of  non-linear  problems  can  be  tackled  in  various  different 
manners.  The  use  of  a  Newton  iteration  approach  involves  linearizing  the  problem  about  the  current  state 
and  inverting  the  global  Jacobian  using  an  appropriate  linear  solver.  Although  Newton  iteration  strategies 
are  widespread,  and  are  amenable  to  modular  software  development  (for  example  in  the  selection  of  existing 
linear  solvers),  for  very  large  systems  of  equations,  the  formation  of  the  Jacobian  matrix  can  be  tedious  and 
memory  intensive.  An  alternate  approach  which  avoids  the  construction  of  the  Jacobian  matrix  consists  of 
using  a  non-linear  iterative  solver  applied  directly  to  the  non-linear  equation  set.  This  approach  has  often 
been  described  as  “placing  the  non-linearity  on  the  inside  of  the  iterative  procedure”  as  opposed  to  “placing 
the  non-linearity  on  the  outside  of  the  iterative  scheme” ,  for  the  former  linearized  approach.  Intuitively,  the 
non-linear  solver  approach  results  in  the  non-linear  residuals  being  updated  frequently  (at  every  iteration), 
whereas  the  linear  solver  approach  results  in  less  frequent  non-linear  residual  evaluations. 

Multigrid  methods  can  be  applied  to  non-linear  problems  either  as  a  linear  solver,  operating  on  a  lin¬ 
earization  of  the  governing  equations,  or  as  a  non-linear  multigrid  formulation  known  as  the  full  approxima¬ 
tion  storage  (FAS)  scheme  [1].  The  fact  that  the  FAS  multigrid  scheme  can  be  applied  directly  to  non-linear 
problems  remains  a  great  advantage  of  multigrid  methods,  which  can  be  used  to  solve  such  problems  in  a 
“matrix  free”  manner.  Within  the  context  of  a  multigrid  method,  a  suitable  error  smoother  must  be  chosen 
to  reduce  high-frequency  errors  on  individual  grid  levels  of  the  multigrid  sequence.  In  the  linear  multigrid 
case,  each  individual  grid  problem  will  be  linear,  while  in  the  FAS  case,  each  grid  problem  will  itself  consist 
of  a  non-linear  problem.  In  this  latter  case,  the  non-linear  smoother  can  either  be  constructed  as  a  non-linear 
iteration  procedure,  or  as  a  linear  solver  operating  on  the  linearization  of  the  local  non-linear  problem. 

The  objective  of  this  paper  is  to  investigate  the  efiectiveness  of  these  various  multigrid  approaches  to 
solving  non-linear  problems.  In  order  to  provide  a  meaningful  comparison,  similar  discretization  and  solution 

*ICASE,  Mail  Stop  132C,  NASA  Langley  Research  Center,  Hampton,  VA  23681-2199,  U.S.A.,  diinitri@icase.edu.  This 
research  was  partially  supported  by  the  National  Aeronautics  and  Space  Administration  under  NASA  Contract  No.  NASl-97046 
while  the  author  was  in  residence  at  ICASE,  NASA  Langley  Research  Center,  Hampton,  VA  23681-2199.  Partial  support  was 
also  provided  under  U.S.  Department  of  Energy  subcontract  B347882  from  Lawrence  Livermore  National  Laboratory. 


1 


strategies  must  be  used  in  the  different  variations  of  the  multigrid  algorithms. 

Our  test  problem  consists  of  a  set  of  non-linear  diffusion  equations,  known  as  the  radiation-diffusion 
equations.  These  equations  govern  the  evolution  of  photon  radiation  in  an  optically  thick  medium,  and 
can  be  derived  from  first  principles  by  integrating  over  all  energy  frequencies,  under  the  assumptions  of 
isotropy  and  small  mean-free  photon  paths  [9] .  These  equations  are  important  in  the  simulations  of  inertial 
confinement  fusion  and  astrophysical  phenomena. 

From  the  numerical  standpoint,  the  radiation-diffusion  equations  provide  a  suitable  test-bed  for  ex¬ 
amining  the  effectiveness  of  numerical  algorithms  at  handling  non-linearities,  since  these  equations  exhibit 
strong  non-linear  behavior,  while  at  the  same  time  they  are  of  modest  complexity,  enabling  a  relatively 
straightforward  Jacobian  construction. 

Various  radiation-diffusion  solution  algorithms  incorporating  multigrid  techniques  have  been  demon¬ 
strated  recently,  but  these  generally  have  employed  linear  multigrid  methods  as  preconditioners  within  a 
Newton-Krylov  method  [2,  4,  8].  Furthermore,  these  multigrid  preconditioners  have  generally  been  applied 
in  the  context  of  an  operator  split  algorithm. 

In  this  work,  we  employ  multigrid  to  the  fully  coupled  system  of  equations.  The  developed  algorithms 
employ  multigrid  directly  as  a  solver  for  the  radiation-diffusion  equations  rather  than  as  a  preconditioner. 
While  the  development  of  preconditioned  Newton-Krylov  solution  techniques  for  these  types  of  problems 
remains  of  interest,  a  more  straight-forward  comparison  between  linear  and  non-linear  multigrid  methods  is 
afforded  when  both  of  these  techniques  are  formulated  as  solvers. 

There  are  various  approaches  to  implementing  multigrid  methods  on  unstructured  meshes.  For  the 
purposes  of  this  comparison,  the  logistics  of  the  multigrid  implementation  are  not  necessarily  the  dominant 
concern,  provided  both  linear  and  non-linear  approaches  are  formulated  in  a  similar  manner.  At  the  outset, 
algebraic  multigrid  methods  cannot  be  considered  herein,  since  these  are  exclusively  formulated  as  linear 
solvers.  Because  our  long  term  objective  is  the  development  of  a  three-dimensional  solver  for  complex 
geometries,  we  avoid  nested  geometric  multigrid  methods,  which  assume  the  existence  of  a  complete  sequence 
of  coarse  and  fine  unstructured  meshes  with  nested  cell  structures,  the  construction  of  which  may  not 
be  feasible  in  the  general  three-dimensional  case.  We  concentrate  instead  on  the  agglomeration  multigrid 
approach,  which  generates  coarse  (non-physical)  levels  automatically  through  a  graph  algorithm  and  makes 
use  of  the  Galerkin  projection  for  constructing  the  coarse  level  equations.  This  approach  can  be  formulated 
in  a  linear  or  non-linear  manner  and  has  been  demonstrated  for  large  three-dimensional  fluid-dynamics 
problems  in  complicated  geometries  [7]. 

The  remainder  of  the  paper  is  organized  as  follows.  In  section  2  we  illustrate  the  correspondence 
between  linear  and  non-linear  multigrid  methods.  The  general  agglomeration  multigrid  approach  as  well  as 
the  three  specific  multigrid  schemes  implemented  for  comparison  are  described  in  section  3.  In  section  4  the 
discretization  of  the  governing  equations  and  sample  test  problem  are  described,  while  section  5  discusses 
our  results,  which  are  also  summarized  in  the  conclusion  in  section  6. 

2.  Linear  and  Non-Linear  MG  Formulations.  The  goal  of  any  multigrid  method  is  to  accelerate 
the  solution  of  a  fine  grid  problem  by  computing  corrections  on  a  coarser  grid  and  then  interpolating  them 
back  to  the  fine  grid  problem.  Although  this  procedure  is  described  in  a  two  grid  context,  it  is  applied 
recursively  on  a  complete  sequence  of  fine  and  coarser  grid  levels.  To  apply  a  linear  multigrid  method  to  a 
non-linear  problem,  a  linearization  must  first  be  performed.  Thus,  if  the  equations  to  be  solved  are  written 
as 
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(2.1) 


r^(Wexact)  —  0 


with  the  current  estimate  w  yielding  the  non-linear  residual  r: 


(2.2)  R(w)  =  r 

the  Newton  linearization  of  this  system  is  taken  as 


(2.3) 


9wh 


Awh  =  — r 


This  represents  a  linear  set  of  equations  in  the  solution  variable  Awh  (the  correction),  to  which  a  linear 
multigrid  (i.e.  MG  correction  scheme)  can  be  applied.  In  this  case,  the  coarse  grid  equation  reads: 


(2.4)  —^Awh  =  -/f  riinear 

aWH 

where  H  and  h  represent  coarse  grid  and  fine  grid  values,  respectively,  and  represents  the  restriction 
operator  which  interpolates  the  fine  grid  residuals  to  the  coarse  grid.  The  residual  of  the  linear  system  on 
the  fine  grid  on  the  fine  grid,  which  is  given  by 


(2.5) 

and  may  be  approximated  as 


I*linear 


gRh 

9wh 


Awh  +  r 


(2.6)  riinear  «  R(w  +  Aw) 

where  R  refers  to  the  non-linear  residual,  as  previously.  The  coarse  grid  corrections  Awh  which  are  obtained 
by  solving  equation  (2.4)  are  initialized  on  the  coarse  grid  as  zero.  After  the  solution  of  equation  (2.4),  these 
corrections  are  prolongated  or  interpolated  back  to  the  fine  grid. 

Alternatively,  a  non-linear  FAS  multigrid  scheme  can  be  used  to  solve  equation  (2.1)  directly  without 
resorting  to  a  linearization.  In  this  case,  the  FAS  coarse  grid  equation  reads: 

(2.7)  Rh('^//)  =  Rh(7^'^/i)  -  7^1* 

where  the  term  on  the  right-hand  side  is  often  referred  to  as  the  defect-correction  [1,  6].  Rh  represents  the 
coarse  grid  discretization  and  and  denote  the  restriction  operators  which  are  now  used  to  interpolate 
residuals  as  well  as  fiow  variables  from  the  fine  grid  to  the  coarse  grids.  In  principal,  different  restriction 
operators  for  residuals  and  variables  may  be  employed.  If  equation  (2.7)  is  re-written  as: 


(2.8) 


Rh('^//)  -Rh(7^'^/i)  =  -Ih^ 
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the  right  hand  sides  of  equations  (2.4)  and  (2.8)  represent  similar  approximations  of  the  restricted  non¬ 
linear  residual,  in  view  of  equation  (2.6)  and  the  fact  that  these  restricted  residuals  in  the  FAS  scheme  are 
always  evaluated  at  the  most  recently  available  fine  grid  updates.  Therefore,  by  equating  the  left  hand  sides 
of  equations  (2.4)  and  (2.8),  the  equivalence  between  the  linear  multigrid  scheme  and  the  non-linear  FAS 
scheme  is  seen  to  be  given  by: 

(2.9)  ^u{wh)  —  ^  Awh 

which  means  that  the  FAS  multigrid  scheme  corresponds  to  an  approximation  to  a  linear  multigrid  scheme, 
where  the  coarse  grid  Jacobians  are  approximated  by  finite  differencing  the  operator.  Therefore,  in  the  limit 
of  asymptotic  convergence,  i.e.  when  Awh  <<  1,  the  two  methods  should  yield  similar  convergence  rates. 

Note  that  the  above  discussion  involves  no  specification  of  the  coarse  grid  operator  and  Jacobian  con¬ 
struction.  Therefore,  a  fair  comparison  of  linear  versus  non-linear  multigrid  methods  should  utilize  a  similar 
construction  for  both  of  these  quantities  in  the  respective  algorithms. 

3.  Multigrid  Algorithms.  The  three  multigrid  variants  implemented  in  this  work  are  based  on  the  ag¬ 
glomeration  multigrid  strategy.  Agglomeration  multigrid  was  originally  developed  for  finite- volume  schemes, 
and  is  based  on  agglomerating  or  fusing  together  neighboring  fine  grid  control- volumes  to  form  larger  coarse 
grid  control  volumes  as  depicted  in  Figure  3.1.  This  approach  has  since  been  generalized  for  arbitrary  dis¬ 
cretizations  following  algebraic  multigrid  principles.  In  fact,  agglomeration  multigrid  can  be  viewed  as  a 
simplification  and  extension  of  algebraic  multigrid  to  non-linear  systems  of  equations.  The  control- volume 
agglomeration  algorithm  can  be  recast  as  a  graph  algorithm,  similar  to  algebraic  multigrid  methods,  where 
the  “seed”  vertex  initiating  an  agglomerated  cell  corresponds  to  a  coarse  grid  point,  and  the  neighboring  ag¬ 
glomerated  points  correspond  to  fine  grid  points,  in  the  algebraic  multigrid  terminology  [10].  While  weighted 
graph  algorithms  can  be  employed  for  agglomeration,  these  weights  cannot  depend  on  solution  values,  as  in 
the  algebraic  multigrid  case,  but  only  on  grid  metrics.  In  this  manner,  the  coarse  grid  levels  are  static  and 
need  only  be  constructed  at  the  beginning  of  the  simulation.  This  avoids  one  of  the  problems  associated 
with  algebraic  multigrid  applied  to  linearizations  arising  from  non-linear  problems,  where  newly  constructed 
coarse  levels  may  be  required  at  each  non-linear  update,  a  task  which  can  be  complicated  in  parallel  com¬ 
puting  environments  [3].  In  the  present  work,  for  simplicity  we  limit  ourselves  to  isotropic  coarsening  using 
an  unweighted  graph,  which  produces  coarse  level  graphs  which  are  maximal  independent  sets  of  the  fine 
grid  graph. 


Fig.  3.1.  Illustration  of  Agglomeration  Multigrid 
Coarse  Level  Construetion 
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As  in  the  algebraic  multigrid  case,  agglomeration  multigrid  employs  a  Galerkin  projection  for  the  construction 
of  the  coarse  grid  equations.  Thus,  the  coarse  grid  operator  is  given  by: 

(3.1)  Rh  =  Ih^hlH 

where  is  the  restriction  operator,  and  is  the  prolongation  operator,  and  both  operators  are  taken  as 
piecewise  constants.  This  simple  construction  applies  equally  to  linear  and  non-linear  operators,  and  reduces 
to  forming  the  coarse  grid  equation  at  an  agglomerated  cell  as  the  sum  of  the  fine  grid  equations  at  each  fine 
grid  cell  contained  in  the  coarse  grid  cell.  The  non-linearities  in  the  operator  are  evaluated  using  solution 
variables  on  the  coarse  grid  interpolated  up  from  the  fine  grid. 

Given  this  multigrid  infrastructure,  three  particular  algorithms  which  differ  mainly  in  the  manner  in 
which  non-linearities  are  handled  are  developed  for  comparison.  The  first  involves  a  standard  non-linear 
FAS  multigrid  algorithm,  and  the  second  involves  a  linear  multigrid  algorithm  applied  to  the  linearization  of 
the  governing  equations.  Finally,  a  hybrid  algorithm  is  also  devised  which  uses  a  non-linear  FAS  multigrid 
outer  cycle,  but  a  linear  solver  on  each  grid  level  as  a  smoother. 

3.1.  FAS  Scheme.  In  the  non-linear  FAS  multigrid  algorithm,  equation  (2.1)  is  solved  directly.  The 
coarse  grid  equations  are  formed  by  Galerkin  projection  (c.f.  equation  (3.1))  and  the  non-linearities  in  the 
coarse  grid  operator  are  evaluated  using  coarse  level  solution  variables  interpolated  up  from  the  fine  grid 
using  the  restriction  operator  (as  per  equation  (2.8)).  On  each  grid  level,  the  discrete  equations  are 
solved  using  a  non-linear  block  Jacobi  iteration  given  as: 


(3.2)  -h  [T)]-^R(w°‘^) 

where  [D]  represents  the  block  diagonal  of  the  Jacobian  matrix.  This  smoother  constitutes  a  non-linear 
solver,  since  the  non-linear  residual  is  updated  at  each  stage,  and  incurs  minimum  memory  overheads  since 
only  the  storage  of  the  block  2X2  matrix  [D]  representing  the  coupling  between  the  two  solution  variables 
at  each  grid  point  is  required. 

3.2.  Linear  Multigrid  Scheme.  The  linear  multigrid  scheme  solves  equation  (2.3)  on  the  fine  grid, 
and  equation  (2.4)  on  the  coarse  levels.  On  the  fine  grid,  the  Jacobian  is  formed  by  explicitly  differenti¬ 
ating  (hand  coding)  the  discrete  operator  Rh-  On  the  coarse  levels,  for  consistency  with  the  FAS  multigrid 
algorithm,  the  Jacobian  is  taken  as  the  explicit  differentiation  of  the  coarse  non-linear  operator  obtained 
by  Galerkin  approximation  (c.f.  equation  (3.1)).  Thus  flow  variables  as  well  as  residuals  are  restricted  to 
the  coarser  grids,  but  the  non-linear  residuals  on  these  coarser  levels  are  not  evaluated,  only  the  Jacobians 
corresponding  to  the  linearization  of  the  non-linear  coarse  level  residuals.  These  coarse  level  Jacobians  are 
evaluated  at  the  beginning  of  the  solution  phase  for  the  non-linear  time-step  problem,  and  are  then  held 
fixed  throughout  the  linear  multigrid  iterations.  Memory  requirements  for  the  linear  multigrid  scheme  are 
increased  over  those  of  the  FAS  scheme  due  to  the  required  storage  of  the  fine  and  coarse  level  Jacobians. 

An  outer  Newton  iteration  is  employed  to  solve  the  complete  non-linear  problem  R(ic)  =  0.  Within  each 
Newton  iteration,  the  linear  system  defined  by  equation  (2.3)  is  solved  by  the  linear  multigrid  algorithm. 
This  provides  a  new  fine  grid  non-linear  correction  Aw  which  is  then  used  to  update  the  non-linear  residual. 
Given  an  accurate  solution  of  the  linear  problem  for  each  outer  iteration,  we  can  expect  the  non-linear 
Newton  scheme  to  converge  very  rapidly  (quadratically) ,  thus  minimizing  the  number  of  non-linear  updates 
required. 
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On  each  grid  level,  the  linear  multigrid  scheme  solves  the  linear  system  using  a  block-Jacobi  smoother. 
If  the  Jacobian  is  divided  up  into  diagonal  and  off-diagonal  block  components,  labeled  as  [D]  and  [O], 
respectively,  the  Jacobi  iteration  can  be  written  as: 


(3.3)  [D]  Awh”+^  =  -r  -  [O]  Awh” 

where  Awh”  represents  corrections  from  the  previous  linear  iteration,  and  Awh”^^  represents  the  new 
linear  corrections  produced  by  the  current  linear  iteration.  At  each  linear  iteration,  the  solution  of  equation 
(3.3)  requires  the  inversion  of  the  block  2X2  matrix  [D]  at  each  grid  point.  The  linear  corrections  Awh  are 
initialized  to  zero  at  the  first  iteration  on  each  grid  level.  Therefore,  this  linear  iteration  strategy  reduces  to 
the  non-linear  Jacobi  scheme  described  above  in  the  event  only  a  single  linear  iteration  is  employed. 

In  contrast  to  the  non-linear  FAS  multigrid  algorithm,  the  residuals,  jacobians  (i.e.  [D]  and  [O]  terms), 
and  the  variables  interpolated  up  to  the  coarse  grids  are  only  evaluated  at  the  start  of  the  non-linear  iteration, 
and  are  held  fixed  during  all  inner  linear  multigrid  cycles  within  a  non-linear  iteration. 

3.3.  Hybrid  Scheme.  For  comparison  purposes,  a  hybrid  scheme  has  also  been  developed  which  is 
based  on  a  non-linear  FAS  multigrid  scheme,  but  solves  the  linearized  equations  on  each  grid  level  using  the 
block  Jacobi  smoother  described  by  equation  (3.3)  instead  of  the  non-linear  iteration  scheme  of  equation 
(3.2). 

While  this  scheme  solves  the  same  form  of  the  equations  on  all  grid  levels  as  the  linear  multigrid  scheme, 
the  non-linear  residuals,  physical  solution  variables,  and  discrete  Jacobians  are  continuously  updated  with 
each  visit  to  a  new  grid  level  in  the  multigrid  scheme,  and  no  outer  Newton  iteration  is  required.  However,  on 
each  grid  level,  multiple  linear  iteration  are  performed  to  solve  the  non-linear  problem  on  that  grid  level,  and 
the  non-linear  residuals  are  thus  only  evaluated  once  on  each  grid  level  for  each  multigrid  cycle,  regardless 
of  the  number  of  linear  iterations  employed.  On  the  other  hand,  this  algorithm  incurs  the  same  memory 
overheads  as  the  linear  multigrid  scheme  due  to  the  required  storage  of  the  Jacobians. 

4.  Problem  Formulation.  The  non-equilibrium  radiation  diffusion  equations  can  be  written  as 


dE 

dt 


V.{DrVE)  =  aa(T^  -  E) 


(4.1) 


— -v.(Avr)  =  -a„(r4-£;) 

with 

—  7^5  Dr{T,E)  =  -  1  I  I  ’  ^t{T)  = 

o(Ta  +  I  I 

Here,  E  represents  the  photon  energy,  T  is  the  material  temperature,  and  k.  is  the  meterial  conductivity. 
In  the  non-equilibrium  case,  the  non-linear  source  terms  on  the  right-hand-side  are  non-zero  and  govern 
the  transfer  of  energy  between  the  radition  field  and  material  temperature.  Additional  non-linearities  are 
generated  by  the  particular  form  of  the  diffusion  coefficients,  which  are  functions  of  the  E  and  T  variables. 
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In  particular,  the  energy  diffusion  coefficient,  Dr{T,E)  contains  the  term  |^|  which  refers  to  the  gradient 
of  E  in  the  direction  normal  to  the  cell  interface  (in  the  direction  of  the  flux).  This  term  constitutes  a  flux 
limiter,  which  is  an  artificial  means  of  ensuring  physically  meaningful  energy  propagation  speeds  (i.e.  no 
larger  than  the  speed  of  light)  [2,  4,  8].  The  atomic  number  2:  is  a  material  coefficient,  and  while  it  may  be 
highly  variable,  it  is  only  a  function  of  position  (i.e.  z  =  f(x,y)  in  two  dimensions). 

Equations  (4.1)  represent  a  system  of  coupled  non-linear  partial-differential  equations  which  must  be 
discretized  in  space  and  time.  Spatial  discretization  on  two-dimensional  triangular  meshes  is  achieved  by  a 
Galerkin  finite-element  procedure,  assuming  linear  variations  of  E  and  T  over  a  triangular  element.  The  non¬ 
linear  diffusion  coefficients  are  evaluated  by  first  computing  an  average  T  and  E  value  along  a  triangle  edge, 
and  then  computing  the  non-linear  diffusion  coefficient  at  the  edge  midpoint  using  these  averaged  values. 
While  a  standard  finite-element  discretization  requires  computing  triangle  averaged  diffusion  coefficients, 
the  edge-based  approach  can  more  easily  be  reproduced  on  coarse  agglomerated  levels  which  do  not  contain 
triangular  structures,  and  simplifies  the  construction  of  exact  discrete  Jacobians.  Both  edge  and  triangle 
based  discretizations  have  been  implemented  and  tested  with  little  observable  difference  in  accuracy.  The 
gradient  of  E  in  the  Dr  diffusion  coefficient  is  also  taken  as  a  one  dimensional  gradient  along  the  direction  of 
the  stencil  edge.  The  source  terms  are  evaluated  using  the  local  vertex  values  of  E  and  T  exclusively,  rather 
than  considering  linear  variations  of  these  variables. 

The  time  derivatives  are  discretized  as  first-order  backwards  differences,  with  lumping  of  the  mass 
matrix,  leading  to  an  implicit  scheme  which  requires  the  solution  of  a  non-linear  problem  at  each  time 
step.  This  approach  is  first-order  accurate  in  time,  and  is  chosen  merely  for  convenience,  since  the  principal 
objective  is  the  study  of  the  solution  of  the  non-linear  system.  More  sophisticated  time  integration  strategies 
involving  higher-order  time  discretizations  and  variable  time-step  sizes,  which  have  been  implemented  by 
other  researches  in  this  field  [2,  5]  are  planned  for  future  work. 

The  Jacobian  of  the  required  linearizations  is  obtained  by  differentiation  (hand  coding)  of  the  discrete 
non-linear  residual.  Because  the  spatial  discretization  involves  a  nearest  neighbor  stencil,  the  Jacobian  can  be 
expressed  on  the  same  graph  as  the  residual  discretization,  which  corresponds  to  the  edges  of  the  triangular 
grid.  The  initial  guess  for  the  solution  of  the  non-linear  problem  at  each  time-step  is  taken  as  the  solution 
obtained  at  the  previous  time-step. 

The  test  case  chosen  for  this  work  is  taken  from  [8]  and  depicted  in  Eigure  4.1.  We  consider  a  unit  square 
domain  of  two  dissimilar  materials,  where  the  outer  region  contains  an  atomic  number  of  2:  =  1  and  the  inner 
regions  (1/3  <  x  <  2/3),(l/3<  y  <  2/3)  contains  an  atomic  number  of  2:  =  10.  The  top  and  bottom  walls 
are  insulated,  and  the  inlet  and  outlet  boundaries  are  specified  using  mixed  (Robin)  boundary  conditions, 
as  shown  in  the  figure.  This  domain  is  discretized  using  a  triangular  grid  containing  7,502  vertices,  shown 
in  Eigure  4.2.  This  grid  conforms  to  the  material  interface  boundaries  in  such  a  way  that  no  triangle  edges 
cross  this  boundary. 

Eigure  4.2  illustrates  a  typical  simulation  for  this  case.  Incoming  radiation  sets  up  a  traveling  thermal 
front  in  the  material,  the  progress  of  which  is  impeded  by  the  region  higher  atomic  number  z.  At  critical 
times  in  the  simulation,  the  diffusion  coefficients  can  vary  by  up  to  six  orders  of  magnitude  near  the  material 
interfaces,  thus  providing  a  challenging  non-linear  behavior  for  the  multigrid  algorithms.  At  each  physical 
time  step,  a  non-linear  problem  must  be  solved.  It  is  the  solution  of  this  transient  non-linear  problem  at 
a  given  time  step  which  forms  the  test  problem  for  the  three  agglomeration  multigrid  algorithms.  Clearly, 
the  size  of  the  physical  time  step  affects  the  stiffness  of  the  non-linear  problem  to  be  solved,  with  smaller 
physical  time-steps  leading  to  more  rapidly  converging  systems.  The  non-dimensional  time-step  chosen  in 
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this  simulation  was  taken  as  0.01.  This  constitutes  a  rather  large  value  compared  to  those  employed  in 
reference  [8]  (usually  of  the  order  of  10“^)  and  may  have  an  adverse  effect  on  overall  temporal  accuracy,  but 
provides  a  more  stringent  test  case  for  the  multigrid  solvers.  Of  the  order  of  1000  time  steps  are  required  to 
propagate  the  thermal  front  from  the  inlet  to  outlet  boundary  in  the  current  simulation. 


Fig.  4.1.  Sample  test  problem  for  non-linear  diffusion  Fig.  4.2.  Illustration  of  unstruetured  grid  for  non¬ 
equations  linear  diffusion  problem:  7,502  vertiees 


Fig.  4.3.  Illustration  of  solution  for  non-linear  diffusion  problem:  Contours  of  T 


5.  Results.  In  the  case  of  the  non-linear  FAS  multigrid  scheme,  the  multigrid  algorithm  is  used  directly 
to  solve  the  non-linear  equations  at  a  given  time  step.  Thus  a  dual-loop  structure  is  used,  containing  an 
outer  loop  over  the  physical  time  steps,  and  an  inner  loop  over  the  number  of  non-linear  multigrid  cycles. 
A  four  level  multigrid  W-cycle  is  employed  by  the  FAS  scheme,  using  6  non-linear  Jacobi  iterations  on  each 
grid  level.  Figure  5.1  shows  the  convergence  achieved  by  the  non-linear  FAS  multigrid  algorithm  at  the  time 
step  corresponding  to  the  center  frame  in  Figure  5.1,  when  the  thermal  front  has  begun  to  encounter  the 
material  interface,  in  the  presence  of  strong  non-linearities.  The  non-linear  residual  is  reduced  by  8  orders 
of  magnitude  over  9  multigrid  cycles  in  this  case,  for  an  average  reduction  rate  of  0.15. 

In  Figure  5.2,  the  convergence  history  for  the  same  case  is  shown  for  the  hybrid  FAS  multigrid  scheme,  using 
6  linear  iterations  per  grid  level  with  a  4  level  W-cycle.  The  convergence  of  this  case  is  almost  identical  to 
that  of  the  previous  case. 

Figures  5.3  and  5.4  depict  the  convergence  of  the  outer  non-linear  Newton  iteration  and  the  inner  linear 
multigrid  iteration.  The  Newton  iteration  is  seen  to  converge  quadratically  for  this  case,  as  expected,  since 
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the  Jacobian  is  a  consistent  differentiation  of  the  discrete  operator.  The  linear  multigrid  iterations  employ 
a  4  level  W-cycle  with  3  linear  Jacobi  iterations  per  grid  sweep  and  3  multigrid  W-cycles  per  non-linear 
iteration.  The  figure  illustrates  the  slight  increase  in  the  linear  residual  at  the  beginning  of  each  non-linear 
iteration,  each  time  the  linear  system  is  modified  by  the  evolution  of  the  non-linearities.  The  large  spike  at 
the  13th  linear  iteration  corresponds  to  the  beginning  of  a  new  physical  time-step  and  thus  new  non-linear 
problem. 


Fig.  5.1.  Convergence  Rate  for  FAS  Multi¬ 
grid  Scheme 


Fig.  5.2.  Convergence  Rate  for  Hybrid  Scheme 


Fig.  5.3.  Convergence  Rate  of  Outer  New¬ 
ton  Reration  for  Linear  Multigrid  Scheme 


Fig.  5.4.  Convergence  Rate  of  Inner  Linear 
Iteration  for  Linear  Multigrid  Scheme 


In  Figure  5.5  the  convergence  of  the  non-linear  residual  for  all  three  methods  is  plotted  in  terms  of  grid 
sweeps,  where  a  grid  sweep  is  defined  as  a  linear  or  non-linear  Jacobi  iteration  on  the  fine  grid.  This  plot 
indicates  that  all  methods  converge  to  approximately  the  same  level  in  40  to  50  grid  sweeps.  However,  when 
these  results  are  plotted  in  terms  of  cpu  time.  Figure  5.6  reveals  the  higher  efficiency  of  the  linear  multigrid 
solver  which  is  almost  an  order  of  magnitude  faster  than  the  non-linear  FAS  multigrid  solver,  while  the 
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hybrid  scheme  lies  in  between  these  two  extremes.  This  is  principally  due  to  the  lower  number  of  non-linear 
residual  evaluations  required  by  the  linear  and  hybrid  schemes. 


Fig.  5.5.  Comparison  of  Convergence  of 
Non-Linear  Residual  in  Terms  of  Number  of 
Grid  Sweeps  for  Three  Multigrid  Schemes 


Fig.  5.6.  Comparison  of  Convergence  of 
Non-Linear  Residual  in  Terms  of  CPU  Time  for 
Three  Multigrid  Schemes 


6.  Conclusions.  The  relative  performance  of  the  three  multigrid  algorithms  described  above  for  a 
particular  time  step  is  indicative  of  the  behavior  observed  at  all  time  steps  throughout  the  entire  simulation. 
For  cases  where  convergence  to  low  residual  tolerances  are  required,  all  schemes  can  be  expected  to  behave 
similarly  on  a  grid  sweep  basis,  since  the  non-linearities  are  essentially  frozen  in  the  asymptotic  convergence 
regime,  but  the  linear  multigrid  method  becomes  more  efficient  simply  due  to  the  smaller  number  of  non-linear 
function  evaluations  required,  a  conclusion  which  is  further  validated  by  the  intermediate  performance  of  the 
hybrid  scheme.  On  the  other  hand,  in  the  initial  phases  of  convergence,  the  FAS  non-linear  multigrid  approach 
can  achieve  faster  convergence  (particularly  on  a  grid  sweep  basis)  since  the  important  non-linearities  are 
advanced  more  rapidly  as  opposed  to  the  linear  method  which  tends  to  “oversolve”  the  linear  system  in  these 
initial  phases  of  convergence.  Another  advantage  of  the  FAS  scheme  is  its  “matrix-free”  formulation  which 
avoids  the  formation  and  storage  of  potentially  large  Jacobian  matrices.  Note  that  this  advantage  is  lost  in 
the  hybrid  scheme,  which  requires  the  Jacobian  matrix  on  each  grid  level,  although  not  simultaneously  on 
all  levels. 

For  this  particular  problem,  the  linear  multigrid  solver  provides  the  overall  most  efficient  solution  strat¬ 
egy.  In  general,  the  relative  performance  of  linear  versus  non-linear  multigrid  methods  depends  on  a  tradeoff 
between  the  expense  of  the  non-linear  function  evaluations  versus  the  size  and  complexity  of  the  Jacobian 
matrices  required  by  the  linearization.  For  cases  where  exact  Jacobians  cannot  be  easily  constructed,  (for 
example  when  larger  discretization  stencils  are  employed),  the  FAS  multigrid  approach  may  be  the  only 
feasible  multigrid  approach  for  solving  the  system  in  a  fully  coupled  manner. 
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